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Abstract. The evolution of an initial perturbation in Vlasov 
plasma is studied in the intrinsically nonlinear long-time limit 
dominated by the effects of particle trapping. After the pos- 
sible transient linear exponential Landau damping, the evo- 
lution enters into a universal regime with an algebraically 
damped electric field, E oc 1/t. The trick used for the Vlasov 
equation is also applied to the two-dimensional (2D) Euler 
equation. It is shown that the stream function perturbation 
to a stable shear flow decays as t~ 5 ^ 2 in the long-time limit. 
These results imply a strong non-ergodicity of the fluid ele- 
ment motion, which invalidates Gibbs-ensemble-based statis- 
tical theories of Vlasov and 2D fluid turbulence. 
PACS 52.35.Ra, 47.27.-i, 47.15.Ki. 

Part of the challenge facing the theory of turbulence 
is that it is extremely difficult to make exact statements 
about the long-time behavior of a nonintegrable system 
that go beyond the mere consequences of applicable con- 
servation laws. For chaotic systems with a few degrees 
of freedom, there are a few results like this, including 
the little-known Sundman's theorems for the three-body 
problem cf. ||, pp. 49-68] and the famous Kolmogorov- 
Arnold-Moser theory [||. Here an attempt is made to 
draw certain long-time conclusions about the nonlinear 
evolution in a Vlasov plasma and in a 2D ideal fluid. We 
study the dynamics of the relaxation of a generic initial 
perturbation in these systems and derive algebraic damp- 
ing laws for the perturbation. As in the above finite- 
dimensional examples, our continuous findings imply the 
lack of ergodicity, with grave implications for several sta- 
tistical theories of turbulence. 

We start with the Vlasov-Poisson system for the elec- 
tron distribution function f(x,v,t) = fo(v) + f(x,v,t) 
and the electric field E(x,t) — —d x <p(x,t), 



(d t + vd x +Ed v )f = 0, d x E 



fdv-l, (1) 



describing nonlinear plasma waves on a uniform ion back- 
ground. In Eq. (0), the time t is normalized to the in- 
verse plasma frequency u>~^ , and x is measured in Debye 
lengths rr) = v e /ujp e , where v e is the electron thermal 
velocity, the unit for v. The problem has two basic di- 
mensionless parameters, the nonlinearity e ~ f / fo and 
the wavenumber kro of the initial perturbation. 

The original solution of the initial- value problem for 
Vlasov plasma by Landau [gj is strictly linear, meaning 
that e is the smallest parameter of the problem. We 
will not assume either of the parameters e or k small or 
large; instead, the largest, or the only large, parameter 
in our treatment will be time. The long-time limit is 
intrinsically nonlinear, because the linearization of the 



Vlasov-Poisson system fails for t larger than the parti- 
cle bounce time Tb ~ e" 1 / 2 This happens because 
the fluctuations of the distribution function do not de- 
cay, but rather develop free-streaming-type small scales, 
f(x,v,t) ~ f(x — vt,v,0) ~ e, and the nonlinearity, 
d v f I fo( v ) ~ increases secularly with time. 

The previous analytical work on the nonlinear Vlasov 
plasma includes the exact special stationary solutions of 
Bernstein, Greene, and Kruskal (BGK) || and the non- 
stationary theory of O'Neil for t\, < t <C and k <^.l 
Q. O'Neil showed that, due to trapping and phase mix- 
ing, the damping rate of the wave, ~/(t) = 4>/4>, starts 
oscillating about zero with the time scale Tf, and a de- 
creasing amplitude. The currently prevailing conjecture 
is that nonlinear plasma waves, after several such oscil- 
lations, settle to a stationary stable BGK wave. This 
conclusion appears to be backed by numerical simulation 
H , although numerical evidence should not be considered 
conclusive for the long-time limit. More importantly, the 
stability of the nonlinear BGK waves remains an out- 
standing issue. This author is not aware of any single 
example of a stable BGK wave; moreover, all analytically 
written BGK waves appear linearly unstable ||, and the 
only known nonlinear stability criterion 0, dfo/dH < 0, 
where H(x, v) = <j)(x) + v 2 /2 is the particle energy, holds 
for no periodic BGK wave ||, P. 85]. This suggests that 
the Landau damping will not be arrested by nonlinearity; 
however, the nature of the damping will be modified for 
large t > e _1 . 

Our logic is as follows. We assume that the electric 
field decays with time: E(x, t) — > 0, as t — > oo. Then this 
assumption is shown to be self-consistent by calculating 
the actual damping rate, E oc t , instead of the linear 
exponential damping. 

Assume a periodic boundary condition in x with the 
period L, and expand the electric field £ in a Fourier 
series. Then, for k ^ 0, the second Eq. (Ill) yields: 



ikE k (t) = {2WL)- 1 / dv dx f(x,v,t)e- lkx 

J -oo JO 

= (2itL)- 1 / fi(a,b)e- ikx ( a ' b >V dadb. (2) 



In Eq. (||), the variables of integration were changed to 
the Lagrangian variables a and 6, the initial position and 
the velocity of a particle. According to the Liouville theo- 
rem, the Jacobian of this transformation is unity, and the 
distribution function is constant along the particle orbit 
thus reducing to its initial value /i(a, b) = f(a, b, 0). 
Equation (|^) expresses the electric field in terms of the 
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particle orbit x(a,b,t) defined by x = E(x,t) and the 
given initial condition, a problem as difficult as the origi- 
nal Eq. (]l|). However, the integral representation of E in 
terms of the orbit is very useful for studying the long-time 
asymptotic, when the electric field is presumably small, 
and the orbit becomes a motion with a constant veloc- 
ity, x(a,b,t) = U(a,b)t (plus lower-order terms). The 
resulting integral of an oscillatory function, 

E k {t) oc J f t (a,b)e~ lktu{a ' b) dadb, t -> oo, (3) 

for smooth /j, will generally have only two kinds of 
asymptotics. If the gradient of U(a, b) is nowhere zero 
(as, for example, in the linear theory, where U ~ 6), 
then the integral (0) is exponentially small at large t 
(the Riemann-Lebesgue lemma). If, on the other hand, 
U has a stationary point where d a U = dbU = 0, then 
the 0(i -1 / 2 )-vicinity of this point dominates the integral, 
which scales as E oc £ -1 . Below we show that U (a, b) has 
stationary points in the general case, and therefore E 
decays algebraically. 
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FIG. 1. Contour lines of U(a,b) computed for exponen- 
tially (left) and algebraically (right) decaying potentials. The 
computation was done through the "infinite" time t — 20 
and t = 1000, respectively. The bottom is a zoom of the 
top. The presense of multiple extrema (O-points) and sad- 
dles (X-points) of U is apparent. 

The problem of finding the final velocity U as a func- 
tion of the initial condition, for a particle moving in a 
decaying potential, is very similar to chaotic scattering 
||, and, likewise, due to the transient particle trapping, 
the function U(a,b) is quite complex (Fig. Q). We are 
interested in whether U(a,b) is a monotonic function of 
its arguments. The fact that it is not is most transpar- 



ent from the inspection of the particle bouncing at the 
top (Fig. ||) and at the bottom (Fig. ||) of a decaying po- 
tential profile. If the initial potential amplitude is small, 
the bouncing at the bottom is possible only if <fi decays 
sufficiently slowly, e.g., <f> oc et~ a , < a < 2, in order 
that the bounce time tj, oc be less than t. The 

initial, linear Landau damping is exponential, seemingly 
suggesting no bouncing, hence no stationary points of 
U(a,b) and the persistence of the exponential damping. 
However, a simple perturbation analysis of the particle 
motion near the top of an evolving potential hill shows 
that one can always pick initial conditions such that the 
behavior of Fig. takes place. To some confusion, this 
turns out possible only if the spatial extrema of <p(x,t) 
and d t <fi(x,t) do not coincide; that is, if there is more 
than just one wave, a safely generic situation. (The re- 
sult of the left Fig. ft] is for two potential waves. A similar 
computation for one wave shows a smooth U with no sta- 
tionary points.) 




FIG. 2. Particle bouncing in a decaying potential and its 
signature U(a,b). Near the potential top, an increase in the 
initial velocity b can bring the particle to the decaying po- 
tential barrier earlier, when the barrier was higher, and thus 
turn the particle around: U(ao,bo) > 0, C/(oo,foo + 6b) < 0. 
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FIG. 3. The cross-section of U(a, b) for the algebraically 
decaying potential shown in Fig. Bj. Near the bottom of the 
potential well, the particle makes many bounces before being 
released in an essentially random direction. 

In fact, U(a,b) has an infinite number of station- 



2 



ary points (a J , W ') . Upon expanding the particle or- 
bit near such a point at large t, x(a,b,t) — UH + 

UL(a - aif + Ul{b - Vf + W ab (a - aP)(b - &*)] t/2+ 

V 1 Int + W' J ' + 0{t- r ), Eq. (§) yields the electric field at 
large t in terms of the infinite series, 



-ikiUH+Vi In t+W j ) 



k 2 Lt 



ULUI - (W ab ) 2 



1/2 



o 



kt 2 



(4) 



which could in principle pose problems in terms of diver- 
gencies or cancellations. 
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FIG. 4. The adiabatic invariant before (a) and after (b) the 
separatrix crossing in a decaying potential well. 

The series ([|) turns out to be absolutely (exponen- 
tially in j) convergent, because it is possible to analyze 
the accumulation of the stationary points of U. This is 
due to the adiabaticity of the particle motion at large 
time, when the bounce frequency u b oc cj) 1 / 2 oc t- 1 ' 2 is 
much larger than the potential damping rate 4>/4> oc t . 
As a result, the adiabatic invariant J(a,6), the (x,v)- 
plane area inside a nearly closed trapped particle orbit, 
is conserved, and the corresponding angle variable 9 is 
growing with the bounce frequency: 9 = J* uj b dt oc t 1 / 2 . 
Untrapping occurs when the shrinking separatrix of the 
decaying potential, with the area S oc i -1 / 2 , intersects 
the orbit with the conserved area J (Fig. [f]). For a small 
J oc (a — ao) 2 + b 2 , the crossing time t* oc J~ 2 and the 
angle 9* oc J -1 . Following a small change during the sep- 
aratrix crossing p0[ , the adiabatic invariant of the pass- 
ing particle (now defined as twice the phase-space area) is 
conserved again and defines the final velocity \U(a, b)\ ~ 
J(a,b)/(2L). The sign of U, roughly sign(sin#*), de- 
pends on whether the crossing happens in the upper or 
in the lower half-plane of Fig. 0. The width of the steps 



of U is still finite, 59* oc e 



-UJ b t 



oc 



e -l/J(a,b) it ig deter _ 



mined by the exponentially narrow near-separatrix layer, 
where the bounce period 2ir/uj b diverges logarithmically, 
and the adiabaticity does not hold. Thus we obtain the 
approximate analytical expression for the final velocity: 

U(a, b) ~ J (a, b)/(2L) tanh \ e 1/J ^ b) sin J' 1 (a, 6)1 . (5) 

Near the bottom of the well, a = ao, the behavior of 
Eq. (U) is consistent with the numerical result in Fig. [jj. 



Equation (||) also implies the exponentially growing cur- 
vature U bb oc e J near the steps as one moves to the accu- 
mulation point of the U extrema, hence the exponential 
convergence of the series (^) . 

Thus the long-time behavior of the electric field (Q) 
is dominated by a few "strongest" stationary points of 
U(a,b). In addition to the algebraic damping rate, we 
infer as a by-product the spectrum E k oc fc~ 2 , fc <C t, 
implying the development of steps in the electron density 
perturbation d x E. 

We now turn to the different problem of the relaxation 
in 2D ideal inviscid incompressible fluid with the velocity 
v = V^>(a:, y,t) x z described by the Euler equation, 



(dt + v ■ V) 



0. 



(0) 



As in the case of Vlasov plasma, we are interested in the 
long-time relaxation of an initial perturbation ?p(x,y,t) 
imposed on a stable shear flow ifj (x) . The deep analogy 
of this problem with the Landau damping in plasmas has 
been noted ||ll| , [12[ . We will assume a periodic boundary 
condition in y [along the shear flow vq(x) = — ip' (x)]. In 
linear theory, the perturbation of the stream function ip 
is known to be damped, because the reconstruction of 
-0 from the conserved vorticity lo with growing gradients 
involves an integration, 

i>{x, y,t) = / G(x, x', y - y') uj(x', y', t) dx' dy' , (7) 



where G is the boundary-condition-dependent Green's 
function with a discontinuous derivative at (x, y) = 
(x',y'). If the flow is unbounded in the x direction, for 
example, Gk{x,x') = e~\ k ^ x ~ x ^. Unlike plasma waves, 
the damping law of tp is algebraic already in linear the- 
ory: ip oc t~ 2 for monotonic Vg (x) 1 1 3 l^Jl^l and oc t~ x l 2 
for vq(x) with an extremum J15| . 

Similarly to the Vlasov case, linear approximation in 
the Euler equation breaks for large t raising the question 
of the long-time asymptotic. To this end, we use the 
same trick as for the Vlasov equation. Upon applying 
the Fourier transform in y to Eq. (Q) and changing the 
integration variables (x',y') to the Lagrangian variables 
(a, 6), we obtain: 

<Afc= [ G k [x,X(a,b,t)}e- lkY{aAt) uj l (a,b)dadb, (8) 



where uii(a,b) is the total initial vorticity, and (X, Y) 
is the orbit of a fluid element with the initial position 
(a, 6). Consider the case of a smooth and monotonic 
vq(x). Then, for a very small perturbation, the unper- 
turbed orbit (A, Y) ~ (a, b + v {a) t) yields an oscillatory 
integral in a, which is not exponentially small because 
of the derivative discontinuity in Gfe. Changing the vari- 
able a to the monotonic Vo(a) and integrating by parts 
twice then yields ipk(x, t) oc t~ 2 for k ^ 0, the well-known 
linear result. For Vq(x) with a stationary point, the sin- 
gularity of the Green's function does not matter, and a 
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stationary-phase integration over a yields ip oc t in 
agreement with fl5|] . Based on the ordering of terms in 
En. (|^) for the regime with ip oc t~ 2 , Brunet and Warn 
||15| argued that the nonlinearity remains small and does 
not change the damping rate. Such an analysis appears 
superficial, because the accumulation of small nonlinear 
effect in the Euler equation is secular. Our analysis of 
Eq. (||) goes as follows. The flow disturbance of order 
e makes the orbit essentially depend on both a and 6, 
e.g., Y = b + vo(a) t + e J vi(a 1 b, t, e) dt. Thus the in- 
tegral (||) is also oscillatory in b for t > e -1 . This is 
when the nonlinearity comes into effect. Because of pe- 
riodicity, the phase iktY(a, b) has a stationary point in 
b producing the additional factor of i -1 / 2 in the integral 
asymptotic. Finally, for a smooth stable monotonic shear 
velocity profile, a smooth stream function perturbation 
decays as ip oc i~ 5 / 2 for t e _1 . Similarly, for vq(x) 
with an extremum, tp oc t . 

One of the interesting consequences of the nonlinear 
Landau damping concerns ergodicity, the assumption un- 
derlying the Gibbs-ensemblc theories of turbulence in the 
Vlasov-Poisson system @,|l7| and in 2D fluid In 
these theories, the statistical ensemble includes all possi- 
ble permutations of phase-space (fluid) elements with the 
associated distribution function / (vorticity lo), via either 
combinatorial treatment or path integration for the parti- 
tion function. Such analyses predict specific quantitative 
results, such as the final relaxed state, for an arbitrary 
initial condition. In addition to the difficulties with non- 



Gaussian path integrals 21 1, the Gibbs-ensemble theory 
cannot be true in such a generality because of the non- 
linear Landau damping. For example, if the initial con- 
dition is a slightly perturbed stable shear flow v y (x) , the 
zonal velocity v x (x,y,t) will decay as eU 5 / 2 , and the 
zonal displacement of any fluid element, SX = J v x dt, 
is for ever bounded by a small constant: \SX\ < Ce. 
This purely dynamical fact does not follow from con- 
servation laws alone and implies that the fluid element 
motion is strongly non-ergodic. It follows that the ex- 
isting statistical theories do not work, at least for initial 
conditions close to stable shear flows. The same is true 
of the Vlasov-Poisson system, where the velocity change 
for any particle, 5v = J E(x(t),t) dt, is bounded for in- 
finite time, because E oc t , and the convergence of Sv 
is ensured by another power of t coming from the nearly 
uniform motion in the coordinate x ~ Ut, over which E 
is zero-average. Again, heating up a small group of parti- 
cles to arbitrarily high energies in an evolving collisionlcss 
plasma does not contradict conservation laws; however, 
this turns out exactly prohibited by self-consistent dy- 
namics. 

The shear damping of perturbations is not specific to 
plane parallel flows in the Euler equation; quite similar 
results must hold for circular monopole vortices develop- 
ing in the course of long-time turbulent evolution [ p2|]23"| 
and also in the framework of related 2D geophysical fluid 
equations, where the nonlinear Landau damping is the 



mechanism of the turbulence relaxation toward large- 
scale coherent structures. Finally, it appears that the 
decaying 2D turbulence is more about dynamics (vortex 
merger and the nonlinear damping of vortex perturba- 
tions) than statistics. 
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